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Abstract 

This paper is devoted to the analysis of a mathematical model of blood cells produc- 
tion in the bone marrow (hematopoiesis). The model is a system of two age-structured 
partial differential equations. Integrating these equations over the age, we obtain a sys- 
tem of two nonlinear differential equations with distributed time delay corresponding 
to the cell cycle duration. This system describes the evolution of the total cell popula- 
tions. By constructing a Lyapunov functional, it is shown that the trivial equilibrium 
is globally asymptotically stable if it is the only equilibrium. It is also shown that the 
nontrivial equilibrium, the most biologically meaningful one, can become unstable via a 
Hopf bifurcation. Numerical simulations are carried out to illustrate the analytical re- 
sults. The study maybe helpful in understanding the connection between the relatively 
short cell cycle durations and the relatively long periods of peripheral cell oscillations 
in some periodic hematological diseases. 

Keywords: Blood cells, hematopoiesis, differential equations, distributed delay, asymptotic 
stability, Lyapunov functional, Hopf bifurcation. 

1 Introduction 

Cellular population models have been investigated intensively since the 1960's (see, for ex- 
ample, Trucco [SSI El], Nooney [25 , Rubinow [2H] and Rubinow and Lebowitz [IS]) and still 
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interest a lot of researchers. This interest is greatly motivated, on one hand, by the medi- 
cal applications and, on the other hand, by the biological phenomena (such as oscillations, 
bifurcations, traveling waves or chaos) observed in these models and, generally speaking, in 
the living world (Mackey and Glass [121, Mackey and Milton :20i). 

Hematopoiesis is the process by which primitive stem cells proliferate and differentiate to 
produce mature blood cells. It is driven by highly coordinated patterns of gene expression 
under the influence of growth factors and hormones. The regulation of hematopoiesis is about 
the formation of blood cell elements in the body. White and red blood cells and platelets 
are produced in the bone marrow from where they enter the blood stream. The principal 
factor stimulating red blood cell production is the hormone produced in the kidney, called 
erythropoietin. About 90% of the erythropoietin is secreted by renal tubular epithelial cells 
when blood is unable to deliver sufficient oxygen. When the level of oxygen in the blood 
decreases this leads to a release of a substance, which in turn causes an increase in the 
release of the blood elements from the marrow. There is a feedback from the blood to the 
bone marrow. Abnormalities in the feedback are considered as major suspects in causing 
periodic hematological diseases, such as auto-immune hemolytic anemia (Belair et al. [4j and 
Mahaffy et al. [13]), cyclical neutropenia (Haurie et al. [T3]), chronic myelogenous leukemia 
(Fowler and Mackey [H] and Pujo-Menjouet et al. [H]), etc. 

Cell biologists classified stem cells as proliferating cells and resting cells (also called Gq- 
cells) (see Mackey [IHIIE])- Proliferating cells are committed to undergo mitosis a certain 
time after their entrance into the proliferating phase. Mackey supposed that this time of 
cytokinesis is constant, that is, it is the same for all cells. Most of committed stem cells 
are in the proliferating phase. The Go-phase, whose existence is known due to the works of 
Burns and Tannock [5, is a quiescent stage in the cellular development. However, it is usually 
believed that 95% of pluripotent stem cells are in the resting phase. Resting cells can exit 
randomly to either entry into the proliferating phase or be irremediably lost. Proliferating 
cells can also be lost by apoptosis (programmed cell death). 

The model of Mackey [TB] has been numerically studied by Mackey and Rey [i21 and 
Crabb et al. [S] . Computer simulations showed that there exist strange behaviors of the stem 
cell population, such as oscillations and bifurcations. Recently, Pujo-Menjouet and Mackey 
[?f] proved the existence of a Hopf bifurcation which causes periodic chronic myelogenous 
leukemia and showed the great dependence of the model on the parameters. 

In this paper, based on the model of Mackey ^16j, we propose a more general model 
of hematopoiesis. We take into account the fact that a cell cycle has two phases, that 
is, stem cells in process are either in a resting phase or actively proliferating. However, 
we do not suppose that all cells divide at the same age, because this hypothesis is not 
biologically reasonable. For example, it is believed that pluripotent stem cells divide faster 
than committed stem cells, which are more mature cells. There are strong evidences (see 
Bradford et al. {[j) that indicate that the age of cytokinesis r is distributed on an interval 
[r, r] with r > 0. Hence, we shall assume that r is distributed with a density / supported on 
an interval [t, t], with < r < r < -foo. The resulting model is a system of two differential 
equations with distributed delay. A simpler model, dealing with the pluripotent stem cell 
population behavior, has been studied by Adimy et al. [1]. 

Some results about stability of differential equations with distributed delay can be men- 
tioned. In [B], Boese studied the stability of a differential equation with gamma-distributed 
delay. Gamma distributions have the property to simplify the nature of the delay and this 
situation is close to the one with discrete delay. Anderson [2 [5] showed stability results 
linked to the different moments (especially the expectation and the variance) of the distri- 
bution. Kuang [TS] also obtained general stability results for systems of delay differential 
equations. More recently, sufficient conditions for the stability of delay differential equations 
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with distributed delay have been obtained by Bernard et al. [5j. They used some properties 
of the distribution to prove these results. However, in all these works, the authors focused 
on sufficient conditions for the stability, there is no necessary condition in these studies, and 
these results are not applicable directly to the model considered in this paper. 

This paper is organized as follows: in Section [21 we present the model and establish 
boundedness properties of the solutions. In Section [3J we study the asymptotic stability of 
the equilibria. We give conditions for the trivial equilibrium to be globally asymptotically 
stable in Section 13.11 and investigate the stability of the nontrivial equilibrium in Section 
13.21 In Section [H we show that a local Hopf bifurcation occurs in our model. In Section \E[ 
numerical simulations are performed to demonstrate that our results can be used to explain 
the long period oscillations observed in chronic myelogenous leukemia. 



2 The hematopoiesis process: presentation of the model 

Denote by r(f, a) and p{t, a) the population densities of resting an proliferating cells, respec- 
tively, which have spent a time a > in their phase at time t > 0. Resting cells can either 
be lost randomly at a rate S > 0, which takes into account the cellular differentiation, or 
entry into the proliferating phase at a rate f3. Proliferating cells can be lost by apoptosis 
(a programmed cell death) at a rate 7 > and, at mitosis, cells with age a divide in two 
daughter cells (which immediately enter the Go-phase) with a rate g{a). 

The function g : [0,t) — + satisfies g{a) ~ Q ii a < t with < r < t < +cx). Moreover, 
it is assumed to be piecewisely continuous such that g{a)da = +00. The later assumption 
describes the fact that cells which did not die have to divide before they reach the maximal 
age T. 

The nature of the trigger signal for introduction in the proliferating phase is not clear. 
However, the work of Sachs [33 shows that we can reasonably think that it strongly depends 
on the entire resting cell population, that is f3 — f3{x{t)), with 

r+oo 

x{t) = / r{t,a)da, t > 0. 



The function /? is supposed to be continuous and positive. Furthermore, from a reasonable 
biological point of view, we assume that /3 is decreasing with lim2;_+_|_oo (3{x) = 0. This 
describes the fact that the rate of re-entry into the proliferating compartment is a decreasing 
function of the Go-phase population. 

Usually, it is believed that the function /3 is a monotone decreasing Hill function (see 
Mackey [E]), given by 

with (3o > 0, 6 > and n > 0. Po is the maximal rate of re-entry in the proliferating phase, 
9 is the number of resting cells at which f3 has its maximum rate of change with respect to 
the resting phase population, and n describes the sensitivity of the reintroduction rate with 
changes in the population. 

The above parameters values are usually chosen (see Mackey [H]) to be 

5 = 0.05day"\ 7 = 0.2 day"\ /3o = 1.77 day" ^ and n = 3. (2) 

Although an usual value 9 is 9 — 1.62 x 10^ cells/kg, it can be normalized without loss of 
generality when one makes a qualitative analysis of the population. 
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Then r(i, a) and p{t, a) satisfy the system of partial differential equations 
dr dr 



{5 + I3{x{t)))r, a > 0, i > 0, (3) 

dt ' da 



dt da 

• -{j + g{a))p, 0<a<r, t>0, (4) 



with 

r(0, a) = i^{a), a > 0, p(0, a) = r(a), a £ [0, t]. 

The functions v ~ v(a) and F = r(a) give the population densities of cells which have spent 
a time a in the resting and proliferating phase, respectively, at time t = 0; that is the initial 
populations of cells with age a in each phase. 

The boundary conditions of system ©-([I]), which describe the cellular flux between the 
two phases, are given by 



r(t,0) = 2 y g(r)p(t, r)dr, 

p(t,0) = fi{x{i))x{t). 

Moreover, we suppose that lima^+oo a) = and lima^TP(i, o) = 0. 

Let y(t) denote the total population density of proliferating cells at time t\ then 



y{t) = / p{t,a)da, t > 0. 
Jo 

Thus, integrating ^ and ^ with respect to the age variable, we obtain 

-{5 + P{x{t)))x{t) + 2 g{T)p{t,T)dT, (5) 

--iy{t) + P{x{t))x{t) ~ ^ 9{T)p{t, T)dT. 



dx 



dy 
dt 



(6) 



We define a function G by 



Git, a) 




a-t 



g{s)ds ) , if t < a, 



g{s)ds I , if a < t. 



Set 

f{T) 9{T)exp[ - I g{s)ds], r > 0. 



One can check that / is a density function, supported on [r, t] , and / represents the density 
of division of proliferating cells. In particularly, f[T)dT = 1. 

Using the method of characteristics to determine p{t,a), we deduce, from ©-([I]), that 
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the process of hematopoiesis is described by the following system: 

"^"^ = -{d + P(x{t)))x{t) 



dt 



2e"^* / G{t,T)T{T ~t)dT, 



0<t<T, 



+ < 



2/ e-''^f{T)P{xit~T))x{t-T)dT 

+ 2e-'^*^ G{t,T)T{T -t)dT, T<t<T, 

2 1 e-'^''f{T)l3{x{t~T))x{t-T)dT, T<t, 



< i < T, 



^ = -jy{t) + P{xit))x{t) 



(7) 



-ft 



J G{t,T)T{T -t)dT, 



- < 



e-^^f{T)l3{x{t ~ T))x{t ~ T)dT 

+ e-T*^ G{t,T)r{T -t)dT, T<t<T, 
e-l'' f{T)P{x{t - T))x{t - T)dT, T<t. 



One can give a direct biological explanation of system ^ . 

In the equation for the resting cells x{t), the first term in the right-hand side accounts 
for Go-cell loss due to either the mortality and cellular differentiation (S) or introduction in 
the proliferating phase . The second term represents a cellular gain due to the movement 
of proliferating cells one generation earlier. It requires some explanations. First, we recall 
that all cells divide according to the density /, supported on [r, r]. We shall call, in the 
following, new proliferating cells, the resting cells introduced in the proliferating phase at the 
considered time t. When i < t, no new proliferating cell is mature enough to divide, because 
cells cannot divide before they have spent a time r in the proliferating phase. Therefore, the 
cellular gain can only proceed from cells initially in the proliferating phase. When t £ [z, r], 
the cellular increase is obtained by division of new proliferating cells and by division of the 
initial population. Finally, when t >t, all initial proliferating cells have divided or died, and 
the cellular gain is obtained by division of new proliferating cells introduced one generation 
earlier. The factor 2 always accounts for the division of each cell into two daughter cells at 
mitosis. The term e""*"*, with t e [0,t], describes the attenuation of the population, in the 
proliferating phase, due to apoptosis. 

In the equation for the proliferating cells y{t), the first term in the right-hand side accounts 
for cellular loss by apoptosis and the second term is for cellular entry from the Go-phase. 
The last term accounts for the flux of proliferating cells to the resting compartment. 

We set /i :— v{a)da. Then, initially, the populations in the two phases are given by 



x(0) pi 



and 



T{a)da. 



At this point, one can make a remark. Since resting cells are introduced in the proliferating 
phase with a rate (3, then r(0), which represents the population of cells introduced at time 
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t = in the cycle, must satisfy 

r(o) = p{fi)fi. 

Taking into account the inevitable loss of proliferating cells by apoptosis and by division, we 
suppose that T{a) is given by 

e~''-p{^i)^l, ifae[0,r), 
na) = { r_.^^.\ ...... .._.r_^ (8) 



'^^ exp ^- J g{s)ds^ 0{fJ.)f^, if a e [z, r). 



This simply describes that T satisfies dH) (see Webb [35], page 8). With ^ and integrating 
by parts, the initial conditions of system ([7]) become 

x(0) = M, y(0) - /3(M)Ai r fir) f i^^) dr. (9) 



When 7 = 0, we have 

y(0) = /3(m)m / Tf{T)dT. 



Assume that the function x t-^ xP{x) is Lipschitz continuous. It is immediate to show 
by steps that, for all > 0, the system ([7]) under condition ([9]) has a unique nonnegative 
continuous solution {x{t),y{t)) defined on [0, +oo). 

One can notice that problem ([7]) reduces to a system of two delay differential equations, 
with initial conditions solutions of a system of ordinary differential equations. On [0,t], the 
first equation for x{t) in system ([TJ reduces to the ordinary differential equation 

^^-{5+ (3m)))m + mp)i^ f e-^^f{T)dT, o<t<T, ^^^^ 

^(0) = 

and, on [t, r], the second equation reduces to the following nonautonomous delay differential 
equation 

^ ^ = -(<5 + /3(^(t)))</^(t) + 2/3(M)M^^-'^V(T)dT 

+ 2 [ e-''^f{T)P{^{t-T)Mt^T)dT, te[T,T], 



where (pit) is the unique solution of (llOp for the initial condition /x. 

By the same way, the solution y{t) of the second equation in denoted by ipit), is given 
in terms of the unique solution ip{t) of pop . associated with /i, and the unique solution ip{t) 
of (HI]), for t e [0,r]. 

Then, the system ^ can be written as an autonomous system of delay differential equa- 
tions, for t >T, 

dx r 

— = - (^ + P{x(t)))x{t) + 2 j e-^V(r)/3(x(t - T))x{t - T)dT, (12a) 
^ = -jyit) + P{x{t))x{t) - r e-'^f{T)fi{x{t ~ T))x(t - T)dT, (12b) 
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with, for t e [0, t], 

xit)=ipit), y{t)=i:{t). (13) 
The solutions of (|12b|) are given expUcitly by 

y{t) ^ f{T)(^j^ e-'^^'-'^(3{x{s))x{s) d^dT fori>T. (14) 

One can notice that y[t) does not depend anymore on the initial population T{a) after 
one generation, that is when t > t. This can be explained as follows: Cells initially in 
the proliferating phase have divided or died after one generation; hence, new cells in the 
proliferating phase can only come from resting cells x{t). 

On the other hand, one may have already noticed that the solutions of (|12ap do not 
depend on the solutions of (jl2bp whereas the converse is not true. The expression of y{t) in 
(|14p gives more precise information on the influence of the behavior of x{t) on the stability 
of the solutions y{t). These results are proved in the following lemma. 

Lemma 2.1. Let {x(t),y{t)) be a solution of il^) . //limf^+oo x{t) exists and equals C > 0, 
then _ 

'1 - e^''^\ 

dr, if 7 > 0, 

1/7 = 0. 




, lim yit) ^{ (15) 



If x{t) is P-periodic, then y[t) is also P-periodic. 
Proof. By using (|14p. we obtain that 



y{t)^ I /(t)( / e~'"Pixit-s))x(t-s) ds]dT for i > r. (16) 

Hence, 



^ hni^ y{t) = /3(C)C /(r) £ e"^^ ds^ dr, 



and (|15p follows immediately. 

When x{t) is P-periodic, then using (|16p it is obvious to see that y{t) is also periodic with 
the same period. □ 

Lemma 12.11 shows the influence of (|12a|) on the stability of the entire system, since the 
stability of solutions of (|12ap leads to stability of the solutions of (|12b|) . 

Before studying the stability of (|12ap . we prove a boundedness result for the solutions 
of this equation. The proof is based on the one given by Mackey and Rudnicki [22] for a 
differential equation with a discrete delay. 

Proposition 2.1. Assume that 6 > 0. Then the solutions of \12a]) are bounded. 

Proof. Assume that 5 > Q and 2(/J' e~^'^ f{T)dT)f3{Q) > S. Since (3 is decreasing and 
limj;^-i_oo f3{x) = 0, there exists a unique xq > such that 

e-^-f{T)dT)f]{xo)=S 
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and _ 

e-^^f{T)dT ] (3{x) < S for x > xq. (17) 



If 2{J^ e-^^f{T)dT)(3{0) < S, then HT]) holds with xq = 0. Set 



One can check that 

e-^'''f{T)dT^ mux (^/3{y)y^ <5x for a; > xi. (18) 
Indeed, let y G [0,x). li y < xq, then 

e-^^f{T)dT]f3{y)y < 2^ / e-'^ f{T)dr) = < Sx 



and, if y > xq, then 

^ e-''^f{T)dT]P{y)y<Sy<Sx. 



Hence, (HH]) holds. 

Assume, by contradiction, that limsupj^^^^ x{t) = +oo, where x{t) is a solution of (jl2ap . 
Then, there exists to > t such that 

x{t) < x(io) for t G [fg ^T,to] and a::(to) > 2^1. 
With (HI]), we obtain that 



e-^V(T)/3(x(to - r))x(to - T)dT < Sx{to). 
This yields, with p2a|l . that 

^(io) < -/3(a;(to)))a:(to) < 0, 

which gives a contradiction. Hence, limsupj^^^ a;(t) < +00. □ 

When (5 = 0, the solutions of ()12a|) may not be bounded. We show, in the next proposition, 
that these solutions may explode imder some conditions. However, one can notice, using p6p . 
that the solutions of (|12bp may still be stable in this case. 

Proposition 2.2. Assume that 6 = and 

e-^^f{T)dT > i (19) 

In addition, assume that there exists a; > such that the function x t—^ x[3{x) is decreasing 
for X >x. If fj,>x, then the unique solution x{t) of \12a\] satisfies 

lim x{t) = +00. 

t — > + CxD 
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Proof. One can notice that, if linif^+oc x{t) = C exists, then (|12ap leads to 

(^2 J\-"'' f{T)dT - 1^(0)0 = 0. 

It follows that C = 0. 

Let fj, >x he given. Consider the equation 

= 2|3{p)^l I e-^''f{T)dT - (3{^{t))lp{t) for < t < r (20) 



with (/3(0) = /i. Since the function x ^ xj3{x) is decreasing for x > x, it is immediate that 
every solution ip{t) of (POjl satisfies, for t G [0,r], 



Consider now the problem 

+2 / e-^-f{T)(i{^{t - T))^{t - T)dT, t e [t,t]. 



(21) 



L (^(0 = lf{t), t e [0,t], 

where (pit) is the unique solution of (|20p for the initial condition /x. Then, 



So, there exists £ > such that r + e < r and (p'{t) > for t G [t, r + e). Since ^ < ip{T) < 
ifiir) < (fiiz + e), for t e [r , r + e] , we have 



V'iz + e) > (2/ e-^V(T)dT-lj/3(^(T + e))<^(r + £) 

e-^V(T)dT)/3(^(T + e))^(r + £) 



> 2/ e-^V(T)dr-l /3(^(T + e))^(r + e). 



Condition (jl9p leads to (p'{t + e) > 0. Using a similar argument, we obtain that 

ip'{t)>0 for te[r,r]. 
To conclude, consider the delay differential equation 

x'{t) = 2f e-'f^f{T)P{x{t~T))x{t~T)dT-f]{x{t))x{t) (22) 



with an initial condition given on [r, t] by the solution (p{t) of (PT|) . Using the same reasoning 
as in the previous cases, we obtain that 

x'(r) > 0. 
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We thus deduce that 

x'{t) > for i > 0. 
This completes the proof. □ 

The assumption on the function x i— > x(3{x) in Proposition 12.21 is satisfied for example 
when (3 is given by H]), with n > 1. In this case, we can take x — 9/{n — 1)^/". 

We now turn our attention to the stability of HH). Problem has at most two equi- 
libria. The first one, = (0,0), always exists: it corresponds to the extinction of the 
population. The second one describes the expected equilibrium of the population; it is a 
nontrivial equilibrium E* = {x* , ?/*), where x* is the unique solution of 



f{T)dT-l]l3{x*)^S (23) 



and, from ([7|) and ([9]), 



P(X*)X* r f{T)(^-^^]dT, if7>0, 

- ^ '' ' i^) 

SX* J Tf{T)dT, if 7 = 0. 

Since /3 is a positive decreasing function and limj,^+oo = 0, then the equilibrium E* 
exists if and only if 

0<S<(2f e-'^^/(r)dr-l )/3(0). (25) 



We shall study in Section [3] the stability of the two equilibria Eq and E* . From Lemma 
12.11 we only need to focus on the behavior of the equilibria of p2ap , that is a; = and x = x*, 
to obtain information on the behavior of the entire population. 



3 Asymptotic stability 

We first show that Eq is globally asymptotically stable when it is the only equilibrium, and 
that it becomes unstable when the nontrivial equilibrium E* appears: a transcritical bifur- 
cation occurs then. In a second part, we determine conditions for the nontrivial equilibrium 
E* to be asymptotically stable. 

3.1 Stability of the trivial equilibrium 

In the next theorem, we give a necessary and sufficient condition for the trivial equilibrium of 
p2ap to be globally asymptotically stable using a Lyapunov functional. Concerning definition 
and interest of Lyapunov functionals for delay differential equations, we refer to Hale |13j . 



Theorem 3.1. The trivial equilibrium of the system U^} is globally asymptotically stable if 

{2 j\-^^f{T)dT-i^m<5 (26) 

and unstable if 

S < (2 r e-'^f{T)dT^l]m- (27) 
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Proof. Wc first assume that (|26p holds. Denote by C+ the set of continuous nonnegative 
functions on [0, r] and define the mapping J : C'^ —^ [0, +00) by 

for all (f E C+, where 



We set (see [13]) 



B{x) = / P{s)s ds for aU x>0. 
Jq 



j{^)= hm sup ■^^'''^ '^^'^^ for ^ e C+ , 
t— 0+ t 



where x''^ is the unique solution of (|12ap associated with the initial condition (p G C+, and 
xf{0) =x'fi{t + 0) for 61 e [0, r]. Then, 



j(^) = ^(r)/3(^(r))^(r) 

+ £ e-^V(r) (^(/3(^(T))^(r))' - (/3(^(r - t))(p(t - T))'^dr. 
Using (|12ap . we have 

- (t) = - (<5 + /3(^(t))) ^(t) + 2 ^ e-^V(T)/9(^(r - r)) ^(r - r)dT. 
Therefore, ipSj) becomes 

AV) = -(5 + P{v(T)))p{^(T))ip\T)+ j\~'<^f{T) (p{ip{T))^{T)f 

+2/3(^(t))v.(t)/3(^(t - t))^(t - r) - (/3(^(t - t))^(t - r)) 
= -(<5 + /5(^(r)))/3(^(r))^2(r) + 2(/5(^(r))^(T))' C e-'^^f{T)dT 



(28) 



e-^V(T) /3((^(T))^(r) - /?(vp(t - t))(^(t - r) 



dr. 



Hence, 



J('/^) < -"(<^(r)), 
where the function u is defined, for a; > 0, by 

u{x) — r{x)(3{x)x'^ 

with _ 

r{x) = (5 - ( 2 / e-^^f{T)dT - 1 )/3(x). 



(29) 



Since /3 is decreasing, r is a monotone function. Moreover, (|26p leads to r(0) > 0, and 
lima;_+oo i~{x) = S > 0. Therefore, r is positive on [0, +cx)). 
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Consequently, the function u defined by (|29p is nonnegative on [0, +00) and u{x) = if 
and only if .x = 0. We deduce that every solution of (|12ap . with (p G C"*", tends to zero as t 
tends to +00. 

We suppose now that (l?7|) holds. The linearization of (|12ap around a; = leads to the 
characteristic equation 

Ao(A) ■.= X + 6 + 13(0) - 2/3(0) ^ e-(^+^)V(r)dr = 0. (30) 

We consider Aq as a real function. Since 

= 1 + 2/3(0) r Te-^^+^'^'f{T)dT > 0, 



it follows that Aq is an increasing function. Moreover, ([50]) yields 

lim Ao(A) = —00, lim Ao(A) = +00, 

A — >~OG X — > + QO 

and (j27p implies that 

Ao{0)^S~(2f e-TV(T)dT - lV(0) < 0. 



Hence, Ao(A) has a unique real root which is positive. Consequently, ([50)1 has at least one 
characteristic root with positive real part. Therefore, the equilibrium x = of (|12ap is not 
stable. This completes the proof. □ 

The inequality is satisfied when S 01 j (the mortality rates) is large or when /3(0) 
is small. Biologically, these conditions correspond to a population which cannot survive, 
because the mortality rates are too large or, simply, because not enough cells are introduced 
in the proliferating phase and, then, the population renewal is not supplied. 

Remark 1. One can notice that when 

e-^^J{T)dT < i, 

the trivial equilibrium Eq is the only equilibrium of (|12p and is globally asymptotically stable. 
When _ 

e-^V(r)dr = 1, 

then Eq is globally asymptotically stable ii 6 > 0. When the equality 

e-'^f{T)dT- 1^(0)^6 

holds, one can check that A = is a characteristic root of ([3D]) and all other characteristic 
roots have negative real parts. Hence, we cannot conclude on the stability or instability of 
the trivial equilibrium Eq of (fT2|l without further analysis. However, this is not the subject 
of this paper. 
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3.2 Stability of the nontrivial equilibrium 

We concentrate, in this section, on the equihbrium E* = {x*, y*) defined by Hence, 
throughout this section, we assume that holds, that is 

0<d<(^2 e-'^f{T)dT - 1^/3(0). 

Since S > and f3{0) > 0, imphes, in particularly, that 

e-^-f{T)dT > i (31) 

From Lemma l2.11 we only need to focus on the stability of the nontrivial equilibrium x = x* 
of p2ap . To that aim, we linearize (|12ap around x* . Denote by /3* G M the quantity 

f3* := A fxPix)) _ ^ = Pix*) + x*P'{x*) (32) 
and set u{t) = x{t) — x* . The linearization of (|12ap is given by 

^ = -{5 + f3*Mt) + W* e-^-f{T)u{t - T)dT. 

Then, the characteristic equation is 

A{X) X + S + P* ~2f]* f e-(-^+'^)^/(r)dr = 0. (33) 



One can notice that the function x i— > x[3{x) is usually not monotone. For example, if [3 
is given by ([T]) with n > 1, the function x i— > xj3{x) is increasing for x < 9/in— 1)^/" and 
decreasing for x > 0/{n — 1)^/^. In this case, (3* is nonnegative when x* is close to zero and 
negative when x* is large enough. 

The following theorem deals with the asymptotic stability of E* . 

Theorem 3.2. Assume that i25\] holds. If 

r > 7f ' (34) 

'i-^^J{T)dT + 1 



then E* is locally asymptotically stable. 

Proof. We first prove that the equilibrium x = x* is locally asymptotically stable when 
/3* > 0. We consider the mapping A(A), given by ([SS]) . as a real function of A. Then A (A) is 
continuously differentiable on K and its first derivative is given by 



dA 

rfA 



1 + 2(3* Te-(^+'^)^/(T)dT > 0. (35) 



Hence, A(A) is an increasing function of A satisfying 



lim A(A) = — oo and lim A(A) +oo. 

A — > — oo A — > + CxD 
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Then, there exists a unique Aq S M such that A(Ao) ~ 0. Moreover, since 

AiO) = S- (^2 J\-''^fiT)dT~iy*, 
we deduce, by using dHB) and that 

A(0) ^-(2 f e-^^f{T)dT - l\x*P'{x*) > 0. 



Consequently, Ao < 0. 

Let X = jj, + iuj he a characteristic root of (|33p such that /i > Aq . Considering the real 
part of we obtain that 



^= -(5 + /?*) + 2/3* y e-(^+T)^/(T)cos(wr)dr. (36) 
Using ([55]) . with A = Aq, together with ((5^ . wc then obtain 

M - Ao = 2/3* / e-'^^fir) [e"''^ cos(cjt) - e"^"^] dr. 



However, 

e"'"^cos(wr) - e"^"^ < 

for all T G [r, r]. So we obtain that — Aq < 0, which leads to a contradiction. This implies 
that all characteristic roots of (|33p have negative real part and the equilibrium x = x* of 
(|12ap is locally asymptotically stable. 
Now, assume that /3* < and 

r > 7¥ • (37) 

e-^''f{T)dT + 1 



Let A = yU + icj be a characteristic root of (|33p such that /i > 0. Since 



dr > 0, 



we have 

2(3*1 e"(^+'')^/(r)cos(wr)dr < -2/3* / e"''^/(r)dT 



So, (IMl) and ([STD lead to 

A* < -((5 + /3*) - 2/3* / e-^'^f{T)dT < 0, 



a contradiction. Therefore, /i < 0. 

Suppose now that p3p has a purely imaginary characteristic root iw, with G M. Then, 
dSi) leads to 

e-''^/(r) cos(wr)dr = '^^f . 
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However, 



■"/(r) cos(wr)(ir 



< / e-''^f{T)dT 



and ([37]) yields 

6 + /3* 



2(i* 



Hence, p3p has no purely imaginary root. Consequently, all characteristic roots of (f55| have 
negative real part and the nontrivial equilibrium a; = a;* of (|12ap is locally asymptotically 
stable. 

Finally, assume that 

P* = ■ (38) 

e~^^f{T)dT + 1 



Consider a characteristic root \ — + itu of pa]), which reduces, with ([38]), to 

\-2f3* e-^V(T)(l + e"^^)dT = 0. (39) 
Suppose, by contradiction, that /i > 0. By considering the real part of pop . we have 
H = 2(3* ( e-^^^fir) (l + e"^^ cos(tJT)) dr < 0. 



We obtain a contradiction, therefore ^ < 0. If we suppose now that /i = then we easily 
obtain that 

cos{ujt) = —1 for all t G [r, r], 

which is impossible. It follows that all characteristic roots of p3p have negative real parts 
when p8|) holds and the equilibrium a:: = x* is locally asymptotically stable. 

From Lemma [2.11 we conclude that E* is locally asymptotically stable when ([34]) holds. 

□ 

The asymptotic stability of E* is shown in Fig. [TJ Values of the parameters are given by 
([!]), except n = 2.42, t — and t = 7 days. The function / is defined by 

0, otherwise. 



The Matlab solver for delay differential equations, dde23 [35], is used to obtain Fig. [T] 
as well as illustrations in Section [U and Section [5] 

When ([34]) does not hold, we have necessarily /3* < 0. In this case, we cannot obtain the 
stability of E* for all values of (3* . In fact, in the next section we are going to show that the 
equilibrium E* can be destabilized, in this case, via a Hopf bifurcation. 
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Figure 1: The solutions x{t) (solid curve) and y{t) (dashed curve) of system (fT2| are drawn 
for values of the parameters /3o, S and 7 given by n — 2A2, t — and t — 7 days. 
In this case, the nontrivial equilibrium E* is locally asymptotically stable, even though the 
solutions oscillate transiently. 



4 Hopf bifurcation and periodic solutions 

In this section, we show that the equilibrium x = x* oi p2ap can become unstable when 
does not hold anymore. Throughout this section, we assume that 

T = 



and 



holds, that is 



0<5<\^2j e-'^^f{T)dT - lj/3(0). 

From Proposition [5Tll the solutions of (|12ap are bounded. Consequently, instability in (|12ap 
occurs only via oscillatory solutions. 
We assume that 

P* < S. (41) 

3-''^/(r)dr + 1 



Otherwise, the nontrivial equilibrium x = x* oi ()12ap is locally asymptotically stable (see 
Theorem 

If instability occurs for a particular value (3* < S, a, characteristic root of must 
intersect the imaginary axis. Hence, we look for purely imaginary characteristic roots iu>, 
w e M, of ([55]) . If iu! is a characteristic root of ([55]) . then is a solution of the system 



d + l3*{l~2Ciu)) = 0, 
uj + 2P*S{lu) = 0, 



(42) 
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where 

C{uj):= [ e-''^/(r) cos(cJT)dT and S{lu) -.^ f e-^^ f{T)sin[uJT)dT. 

Jq Jq 

One can notice that = is not a solution of ([1^ . Otherwise, 



5^ (^2j\-^^f{T)dT~?jP*<0, 



which gives a contradiction. Moreover, if is a solution of (|42p . then —iuj is also a charac- 
teristic root. Thus, we only look for positive solutions lu. 

Lemma 4.1. Assume that the function t i— > e^'^'^ fir) is decreasing. Then, for each S such 
that \25\) is satisfied, has at least one solution [f3*,ijJc), with (3* < S and > 0. It 

follows that f33\) has at least one pair of purely imaginary roots ±.iujc for (3* — l3*. Moreover, 
iitOc are simple characteristic roots of i33\) . Consider the branch of characteristic roots 
A(— /3*), such that A(— /?*) = iuc- Then 



dRe{X) 



> if and only if - S[ ] > C'K). (43) 



Proof. First, we show by induction that S{uj) > for uj > 0. It is clear that S{uj) > if 
OJT G (0, tt]. Suppose that lot G (tt, 27r]. Then 



S(lo) = - / e-'if(-]sm(T)dT 
^ Jo Vcj/ 

= - / e-'^-zf-) sin(T)dT + - / e-^^ f(-\sin(T)dT. 
UJ Jq \loJ uj J^ \ujJ 

Since / is supported on the interval [0,t], it follows that 

£e-^/(9sin(r)dr = 0. 



So, we obtain 



Since the function t t-^ e~^'^ f{T) is decreasing, we finally get S{uj) > 0. Using a similar 
argument for ujt G {kir, [k + l)7r], with fc G N, fc > 2, we deduce that S{uj) > for all a; > 0. 
Consider the equation 

, , uj(l-2C(uj)) , , 



The function g is continuous with 

Te-''^f{T)dT 



hm = < (45) 
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because ([25]) leads to 1 — 2C(0) < 0. Moreover, the Riemann-Lebesgue's lemma implies that 



lim C{uj) = lim S{lu) = 0. 

uj — >+oo ■ ■ ■ ' 



This yields 



lim .g(cij) = +0O. 



UJ — ^ + oo 



We conclude that there exists a solution ujc > oi 
we obtain 1 - 2C{ujc) > 0. Set 



Since |C(wc)| < C(0), it follows that 



1 - 2C(wc) 



2C(0) + 1 



< 



Since S{uJc) > and g{ujc) — S > 0, 

(46) 



One can check that {/3*,ujc) is a solution of (|^^ . It follows that ztiujc are characteristic roots 
of® for /3* 

Define a branch of characteristic roots A(— /?*) of such that A(— /3*) = io^c- We use 
the parameter —(3* because (3* < 6 < 0. 
Using we obtain 



1 + 2/?* / re 



-(A+7)t 



/(r)dr 



dX 



1 - 2 



-(A+7)t 



/(r)dr. 



d(-/3*) 

If we assume, by contradiction, that iuJc is not a simple root of then ([i7|) leads to 

CiiOc) = ^ and ^(wc) = 0. 

Since S{ujc) > 0, we obtain a contradiction. Thus, itOc is a simple root of 
Moreover, using (I47p . we have 



(47) 



1 + 2/3* / re 



^(A+7)r 



/(r)dr 



1-2 / e-(^+^)V(r)dr 
Jo 

Since A is a characteristic root of we also have 

A + 5 



1 - 2 



-(A+7)t 



/(r)(ir 



/3* 



So, we deduce 



dX 



1 + 2/3* / re 



-(A+7)t 



/(r)dr 



A + (5 
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Then, 



sign 



f dRejX) \ 
\di-p*)j 



sign< Re 



f3'=t3' 



dX 



= sign<^Re 
= signl - (3, 



1 + 2/3* f Te-(^+'f^^f{T)dT 

Jo 

X + 6 

Sjl + 2p*S'iu;,)) + 2(3*uJcC'{lUc) 

(52 +tj2 



= sign^S{l + 2/3*S'{uJc)) + 2(3lujcC' (uJc) 
From and the fact that 1 — 2C{ijJc) > 0, this leads to 

= signl 1 - 2C{uJc) - 2SS'{ujc) - 2lUcC'{lUc) 



(dRe(X)] 



Sign 



l^2u;c(^ 



S{uJc)\' 



S{0Jc)\' 



= sign<^ — C'{uJc) — 5 

This concludes the proof. 

Remark 2. Consider the function g defined by (|44p and denote by a the quantity 



a:= \ 2 e-^V(T)dT-lj/3(0). 

Define the sets 

f2 := {w > 0; < .g(w) < a and c/'(tj) = 0} and A := 
One can notice that A is finite (or empty). If 5 G (0, a) \ A, then 

dRe{X) 



d{-f3*) 



^0. 



Indeed, we have 



S{lo) 



w > 0. 



□ 



Since (5 ^ A, we have g'{LOc) 7^ 0. Moreover, g{oJc) = 5- Thus 

We conclude by using ([43]) . 

Lemma HTI together with Remark 2, allows us to state and prove the following theorem. 
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Theorem 4.1. Assume that the function t i~-> e^'^'^/(r) is decreasing. Then, for each 6 ^ A 
satisfying i25\) . there exists (3* < 5 such that the equilibrium x = x* is locally asymptotically 
stable when [3* < (3* < 5 and a Hopf bifurcation occurs at x = x* when (3* = (3*. 

Proof. First, recall that x = x* is locally asymptotically stable when /3* ^ S (see Theorem 
I3.2p . We recall that, from the properties of the function g, equation (j44|l has a finite number 
of solutions (see Lemma HTTj) . We set 



where u* is the smaller positive real such that 

C{lu*) = iam{C{Lu)] w is a solution of (gH)}. 

Then, /3* is the maximum value of /3* (as defined in Lemma l4.ip which gives a solution of 
(jUJ. From LemmaHHl (|33p has no purely imaginary roots while (3* < (3* < 5. Consequently, 
Rouche's Theorem [TOl p. 248] leads to the local asymptotic stability of a; = a;*. 

When (3* = f3*, (|33|) has a pair of purely imaginary roots UlUc, uJc > (see Lemma [4. ip . 
Moreover, since S ^ A, Remark 2 implies that 



dReiX) 



Assume, by contradiction, that 



dRe{X) 



^0. 

=/5* 



< 



d{-(3*) 

for j3* > (3*, (3* close to (3*. Then there exists a characteristic root \{—(3*) such that 
ReX{—(3*) > 0. This contradicts the fact that a; = a;* is locally asymptotically stable when 
(3* > (3*. Thus, we obtain 

dRe{X) 

This implies the existence of a Hopf bifurcation a,t x = x* for (3* — (3*. □ 



> 0. 



With the values of 5, 7 and (3o given by and t — 1 days, has periodic solutions 
for /3* = —0.3881 with a period about 33 days. This value of /?* corresponds to n = 2.53 (see 
Fig. ^andEJ. The function / is given by (liU]) . 

The bifurcation parameter was chosen to be /3* in this study, and the values of /?* depend 
strongly on the sensitivity n of the function /3(x), since all other parameters are fixed by ()25p . 
In this model, the sensitivity n plays a crucial role in the appearance of periodic solutions. 
Pujo-Menjouet and Mackey ^\ already noticed the influence of this parameter on system 
p2p when the delay is constant (or equivalently, when / is a Dirac measure) . The sensitivity 
n describes the way the rate of introduction in the proliferating phase reacts to changes in 
the resting phase population produced by external stimuli: a release of erythropoietin, for 
example, or the action of some growth factors. 

Of course, the influence of other parameters (like mortality rates b and 7, or the minimum 
and maximum delays t and r) on the appearance of periodic solutions could be studied. 
However, since periodic hematological diseases — defined and described in Section [5] — 
are supposed to be due to hormonal control destabilization (see .11, ). then the parameter 
n, among other parameters, seems to be appropriate to identify causes leading to periodic 
solutions in ([T^ . 
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Figure 2: The solutions of system p^ . x{t) (solid curve) and y{t) (dashed curve), are drawn 
when the Hopf bifurcation occurs. This corresponds to n = 2.53 with the other parameters 
given by ^ and t = 7 days. Periodic solutions appear with period of the oscillations about 
33 days. 




Figure 3: For the values used in Figl2 the solutions are shown in the (x,y)-plane: the 
trajectories reach a limit cycle, surrounding the equilibrium. 
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5 Discussion 

Among the wide range of diseases affecting blood cells, periodic hematological diseases (Hau- 
rie et al. ) are of main importance because of their intrinsic natme. These diseases are 
characterized by significant oscillations in the number of circulating cells, with periods rang- 
ing from weeks (19-21 days for cyclical neutropenia JA\) to months (30 to 100 days for chronic 
myelogenous leukemia [14]) and amplitudes varying from normal to low levels or normal to 
high levels, depending on the cells types [Tl]. Because of their dynamic character, periodic 
hematological diseases offer an opportunity to understand some of the regulating processes 
involved in the production of hematopoietic cells, which arc still not well-understood. 

Some periodic hematological diseases involve only one type of blood cells, for example, 
red blood cells in periodic auto-immune hemolytic anemia (Belair et al. [1]) or platelets in 
cyclical thrombocytopenia (Santillan et al. [31 ). In these cases, periods of the oscillations 
are usually between two and four times the bone marrow production delay. However, other 
periodic hematological diseases, such as cyclical neutropenia (Haurie et al. [H) or chronic 
myelogenous leukemia (Fortin and Mackey [H]), show oscillations in all of the circulating 
blood cells, i.e., white cells, red blood cells and platelets. These diseases involve oscillations 
with quite long periods (on the order of weeks to months). A destabilization of the pluripo- 
tential stem cell population (from which all of the mature blood cells types are derived) seems 
to be at the origin of these diseases. 

We focus, in particularly, on chronic myelogenous leukemia (CML), a cancer of the white 
cells, resulting from the malignant transformation of a single pluripotential stem cell in the 
bone marrow (Pujo-Menjouet et al. [26]). As described in Morley et al. [24^, oscillations can 
be observed in patients with CML, with the same period for white cells, red blood cells and 
platelets. This is called periodic chronic myelogenous leukemia (PCML). The period of the 
oscillations in PCML ranges from 30 to 100 days (Haurie et al. [T^, Fortin and Mackey [TT|). 
depending on patients. The difference between these periods and the average pluripotential 
cell cycle duration (between 1 to 4 days, as observed in mouses, see Mackey [18]) is still not 
well-understood . 

Recently, to understand the dynamics of periodic chronic myelogenous leukemia, Pujo- 
Menjouet et al. considered a model for the regulation of stem cell dynamics and investi- 
gated the influence of parameters in this stem cell model on the oscillations period when the 
model becomes unstable and starts to oscillate. In this paper, taking into account the fact 
that a cell cycle has two phases, that is, stem cells in process are either in a resting phase or 
actively proliferating, and assuming that cells divide at different ages, we proposed a system 
of differential equations with distributed delay to model the dynamics of hematopoietic stem 
cells. By constructing a Lyapunov functional, we gave conditions for the trivial equilibrium 
to be globally asymptotically stable. Local stability and Hopf bifurcation of the nontrivial 
equilibrium were studied, the existence of a Hopf bifurcation leading to the appearance of 
periodic solutions in this model, with a period around 30 days at the bifurcation. 

Numerical simulations show that periodic solutions occur after the bifurcation, with peri- 
ods increasing as the bifurcation parameter (the sensitivity n) increases. In Fig. |31 solutions 
oscillate around the equilibrium values with periods around 45 days. Moreover, amplitudes 
of the oscillations range from low values to normal values. The sensitivity is equal to n = 3, 
that is the parameters are given by This corresponds to values given by Mackey [TO] , 
values for which abnormal behavior (periodic) is usually observed in all circulating blood 
cells types. 

When n continues to increase, longer oscillations periods are observed with amplitudes 
varying from low values to high values (see Fig. [5]). This situation characterizes periodic 
chronic myelogenous leukemia, with periods in the order of two months (70 days). 
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Figure 4: Solutions x{t) (solid curve) and y{t) (dashed curve) of system oscillate with 
periods close to 45 days; the parameters are the same as in Fig [21 with n = 3. The amplitudes 
of the oscillations range from low values to normal values. 




Figure 5: Solutions x{t) (solid curve) and y{t) (dashed curve) of system IT^ oscillate with 
periods close to 70 days; the parameters are the same as in FigO with n = 4. The amplitudes 
of the oscillations range from low values to high values. 
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Moreover, the oscillations observed in Fig. 2] and [5] look very much like relaxation oscilla- 
tions. Experimental data from patients with PCML suggest that the shape of oscillations is 
of a relaxation oscillator type [TTJ[T3]. Furthermore, Fowler and Mackey [T^] showed that a 
model for hematopoiesis with a discrete delay may also exhibit relaxation oscillations. There- 
fore, it seems that not only periods and amplitudes of the oscillations correspond to the ones 
observed in PCML, but also the shape of the oscillations. 

Numerical simulations demonstrated that long period oscillations in the circulating cells 
are possible in our model even with short duration cell cycles. Thus, we are able to charac- 
terize some hematological diseases, especially those which exhibit a periodic behavior of all 
the circulating blood cells. 
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